function [ s ] = evaluateQuadraturePoint( p, q )
%EVALUATEQUADRATUREPOINT Summary of this function goes here
%   Detailed explanation goes here

% shiftAlpha = 1e-3;

global shiftAlpha shearModulus poissonRatio;

r = p - q;

rn = max(norm(r), shiftAlpha);

rg = r / rn;

% qnc = [qn(2); -qn(1)];

M = [rg(1)*rg(1), rg(1)*rg(2), rg(1)*rg(3); rg(2)*rg(1), rg(2)*rg(2), rg(2)*rg(3); rg(3)*rg(1), rg(3)*rg(2), rg(3)*rg(3)];

I = diag(ones(3, 1));

s = 1 / (16 * pi * (1 - poissonRatio) * shearModulus * (rn + shiftAlpha)) * ((3 - 4 * poissonRatio) * I + M);

end

